Edge species

Sooooo nothing with BEC zones has been promising so far. Let’s try something else. Let’s switch to VRI data and see if amount of mature/old growth forest is related to amount of squirrel in the diet, since squirrels like old forests.

# Load some libraries.
library(tidyverse)
library(landscapemetrics)
library(raster)
library(sf)
library(AICcmodavg)
library(MuMIn)

# Bring in diet data.
df <- read_csv('../data/interim/camera_corrected.csv', guess_max=7000)
# I had some trouble with readr so I increased the number of rows used to guess.

source('../src/prey_attributes.R')

head(items)

# Calculate proportion of diet made up of squirrel at each site.
sq.mass <- items %>% mutate(mass=as.numeric(mass)) %>% 
  group_by(site) %>% 
  mutate(total=sum(mass)) %>% 
  filter(genus == 'Tamiasciurus') %>% 
  mutate(amount.sq=sum(mass), prop.sq=amount.sq/total) %>% 
  dplyr::select(site, prop.sq) %>% distinct()

sq.mass

That’s remarkably consistent.

The next step is to get the amount of mature forest for each of these sites. Unfortunately, the VRI data for Turbid Creek is useless so that drops me down to five sites. But I can pull in the VRI for the others…

# Import transition zone shapefile.
vri <- st_read('../data/external/VRI_camera-sites_2019.shp')

Reading layer `VRI_camera-sites_2019' from data source `C:\Users\gwync\sfuvault\productivity-occupancy\data\external\VRI_camera-sites_2019.shp' using driver `ESRI Shapefile'
Simple feature collection with 13321 features and 188 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 413289.1 ymin: 5428268 xmax: 645675.7 ymax: 5589928
projected CRS:  WGS 84 / UTM zone 10N

Then I need to rasterize it. Which means first I need to break it down into classes so I can assign a raster value. Unfortunately there’s no single field in the VRI data that works for classification. I can break this into:

• Water (1)
• Land (unvegetated) (2)
• Vegetated (not forest) (3)
• Deciduous (4)
• Coniferous (regen, young, medium, and old) (5, 6, 7, 8)

# Assign to classes.
# Regen/young/mature/old ages come from Zharikov et al. 2007
vri.class <- vri %>% mutate(hab.class=case_when(
  # Non-vegetated and water
  BCLCS_LV_1 == 'N' & BCLCS_LV_2 == 'W' ~ 1, 
  # Non-vegetated and not water
  BCLCS_LV_1 == 'N' & BCLCS_LV_2 != 'W' ~ 2, 
  # Vegetated and not trees
  BCLCS_LV_1 == 'V' & BCLCS_LV_2 == 'N' ~ 3,
  # Vegetated and mixed or deciduous trees
  BCLCS_LV_1 == 'V' & BCLCS_LV_4 %in% c('TB', 'TM') ~ 4,
  # Vegetated and coniferous trees
  BCLCS_LV_4 == 'TC' & PROJ_AGE_1 < 20 ~ 5,
  BCLCS_LV_4 == 'TC' & PROJ_AGE_1 >= 20 & PROJ_AGE_1 < 60 ~ 6,
  BCLCS_LV_4 == 'TC' & PROJ_AGE_1 >= 60 & PROJ_AGE_1 < 140 ~ 7,
  BCLCS_LV_4 == 'TC' & PROJ_AGE_1 >= 140 ~ 8,
  TRUE ~ 0
))

# See how it turned out.
vri.class %>% group_by(hab.class) %>% summarise(n())

Simple feature collection with 9 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 413289.1 ymin: 5428268 xmax: 645675.7 ymax: 5589928
projected CRS:  WGS 84 / UTM zone 10N

A little investigating shows that the 4 polygons classed as 0 are unclassified in the VRI or conifer polygons with unknown ages, so that seems fine.

# Set raster extent based on tz shapefile.
ext <- extent(vri)

# Make an empty raster to populate with values.
r <- raster(ext, res=c(100, 100))

# Populate BEC polygon data onto empty raster grid.
r.vri <- rasterize(vri.class, r, 'hab.class')

# Create labels for raster.
vri.levels <- data.frame(ID=0:8, class.name=c('undefined', 'water', 'land', 'unforested', 'deciduous', 'regen', 'young', 'mature', 'old'))

# Add them to the raster.
levels(r.vri) <- vri.levels

# Save the raster image.
writeRaster(r.vri, '../data/interim/vri_camera-sites_2019.tif', format='GTiff')

If picking up later, load up the raster and keep going.

# Import the raster.
r.vri <- raster('../data/interim/vri_camera-sites_2019.tif')

# Import nests and calculate centroids.
sites <- read_csv('../data/processed/the_big_list_of_nests.csv') %>% 
  group_by(name) %>% 
  mutate_at(c('lat', 'lon'), mean) %>% 
  mutate_at(vars(starts_with('status')), max) %>% 
  mutate_at(c('telemetry', 'cameras', 'remains'), max) %>% 
  dplyr::select(-nest, -NOTES) %>% 
  distinct() %>% 
  filter(cameras > 0)

# Drop TCR, ungroup.
sites <- sites %>% filter(site != 'TCR') %>% 
  ungroup()

# Make site table a spatial object and make it UTMs.
sites.sf <- ungroup(sites) %>% st_as_sf(coords=c('lon', 'lat')) %>%
  st_set_crs('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') %>%
  st_transform("+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs")

# Convert approximate homerange area in ha to radius in m
a.hr.ha <- 3700
r.hr.m <- sqrt(a.hr.ha*10000/pi)

# Calculate area per class per site.
class.area <- sample_lsm(r.vri, y=sites.sf, size=r.hr.m, what='lsm_c_ca', 
           shape='circle') %>% 
  # Amend class number with habitat class name.
  left_join(data.frame(levels(r.vri)), by=c('class'='ID')) %>% 
# Reorganize by site.  
  pivot_wider(id_cols=plot_id, names_from=category, values_from=value, 
              values_fill=list(value=0))

# Amend with site names.
class.area <- dplyr::select(sites, site) %>% rownames_to_column() %>% 
  mutate(rowname=as.integer(rowname)) %>% 
  right_join(class.area, by=c('rowname'='plot_id'))

class.area

Ok that looks good so far. Obviously Ruby Lake is the one with all the water. Otherwise, there seems to be a fair bit of variability regarding the different cover types.

I have two age classes that might be relevant to squirrels: mature and old. So amount of squirrel might be related to mature or old, or it might be related to mature and old. (Or, more likely, it is related to neither.)

# Join VRi and diet data together.
sq.vri <- sq.mass %>% filter(site != 'TCR') %>% left_join(class.area)

# Sum mature and old forest classes
sq.vri <- sq.vri %>% mutate(older=sum(mature, old))

sq.vri

# Plot it.
ggplot(sq.vri, aes(x=older, y=prop.sq, label=site)) +
  geom_point() +
  geom_text(hjust='inward', nudge_y = 0.5) +
  geom_smooth(method='lm', se=FALSE, color='black', linetype='dashed') +
  theme_classic()

Shocking! There’s actually kind of almost a little bit of a pattern. So here are my models:

# Proportion of squirrel biomass is function of amount of mature forest.
sq.x.mature <- lm(prop.sq ~ mature, data=sq.vri)

# Proportion of squirrel biomass is function of amount of old forest.
sq.x.old <- lm(prop.sq ~ old, data=sq.vri)

# Proportion of squirrel biomass if function of total older forest.
sq.x.older <- lm(prop.sq ~ older, data=sq.vri)

# Do the AIC thing.
models <- list(sq.x.mature, sq.x.old, sq.x.older)
modnames <- c('mature', 'old', 'older')

aictab(models, modnames=modnames)

Model selection based on AICc:

       K  AICc Delta_AICc AICcWt Cum.Wt   LL
older  3 17.40       0.00   0.62   0.62 6.30
old    3 19.00       1.61   0.28   0.89 5.50
mature 3 20.93       3.53   0.11   1.00 4.54

# And the R2 thing.
summary(sq.x.mature)$r.squared

[1] 0.1190584

summary(sq.x.old)$r.squared

[1] 0.4006557

summary(sq.x.older)$r.squared

[1] 0.5652464

Which does not look promising. But the other option is to look at core area, assuming squirrels avoid edge. The problem here is that I should combine both old and mature forest into a single category, so that edge between them doesn’t count, and that probably requires remaking my raster.